

% function takes image of vessel centers selected in FIJI and draws
% lines for linescans from it; input is just either the filename of the
% manually selected image or the actual image; outputs are ep1 and ep2
% which are variables where the first column is the x-coordinate and second
% column is the y-coordinate of the first endpoint (ep1) and second
% endpoint (ep2) of the lienscan line
% figtitle is an optional input that will display the name of the .tif file
% 
function [ep1,ep2,endpointIDs] = linescansFromCentersMalabSelect(mov,micronsPerPixel,figtitle)

if nargin < 2 || isempty(micronsPerPixel)
    micronsPerPixel = 200/512;      %this corresponds to a zoom of 5
end

%% select vessel centers:
if size(mov,3) > 30
    im = mov(:,:,30);
else
    im = mov(:,:,1);
end
if nargin < 3 || isempty(figtitle)
    figtitle = 'select vessel';
end
[segs,imsz] = selectVesselCenters(im,figtitle);
im = interpolateLineSegs(segs,imsz);        %note this function will return identities of each line segment; currently these aren't being used but could be in the future

vesselIDs = im;          %this will keep vessel IDs if they're available
im(im~=0) = 1;

ind = find(im == 1);
orientIm = nan(size(im));
for n = 1:length(ind)
    orientIm(ind(n)) = getLocalOrientation(im,ind(n));
end
orientIm = deg2rad(orientIm);
%% downsample the number of points
downsampleDistance = 10;
[~,downsampledInd] = downsampleConnected(im,downsampleDistance);

%% draw linescans for all these points:

maxradius = 40/micronsPerPixel;             %can reduce this number to have a shorter linescan (if your vessels are all narrow)
%maxradius = 10/micronsPerPixel;             %can reduce this number to have a shorter linescan (if your vessels are all narrow)
ep1 = zeros(length(downsampledInd),2);
ep2 = ep1;
[M,N] = ind2sub(size(im),downsampledInd);
endpointIDs = zeros(size(ep1,1),1);             %saves the vessel# from which that endpoint came
for n = 1:length(downsampledInd)
    ep1(n,1) = sin(orientIm(downsampledInd(n)))*maxradius + N(n);
    ep1(n,2) = cos(orientIm(downsampledInd(n)))*maxradius + M(n);
    
    ep2(n,1) = -sin(orientIm(downsampledInd(n)))*maxradius + N(n);
    ep2(n,2) = -cos(orientIm(downsampledInd(n)))*maxradius + M(n);
    endpointIDs(n) = vesselIDs(downsampledInd(n));
end

%% testing block to check that linescans are drawn properly
if 1
    figure
    imagesc(orientIm)
    hold on
    for n = 1:size(ep1,1)
        plot([ep1(n,1),ep2(n,1)],[ep1(n,2),ep2(n,2)],'o-')
    end
    axis equal
end



